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Abstract. The objective of this work is to understand how the characteristics of 
relativistic MHD turbulence may differ from those of nonrelativistic MHD turbulence. 
We accomplish this by studying the ideal invariants in the relativistic case and 
comparing them to what we know of nonrelativistic turbulence. Although much work 
has been done to understand the dynamics of nonrelativistic systems (mostly for ideal 
incompressible fluids), there is minimal literature explicitly describing the dynamics of 
relativistic MHD turbulence using numerical simulations. Many researchers simply 
assume that relativistic turbulence has the same invariants and obeys the same 
dynamics as non-relativistic systems our results show that this assumption may be 
incorrect. 
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1. Introduction 

Many studies in numerical relativity and high-energy astrophysics depend on the 
dynamics of relativistic plasmas. These include phenomenon such as primordial 
turbulence, neutron stars, active galactic nuclei, and accretion disks near black 
holes P m El [HI [121 ESI EHl [30]. Unfortunately, we do not know if the results of these 
studies are accurate because of approximations such as the use of nonrelativistic fluid 
dynamics and the lack of a standard model to describe the dynamics of relativistic 
turbulence. In particular, very little is understood about the turbulent dynamics 
of a relativistic plasma or its effect on the evolution of magnetic fields. This can 
only effectively be studied through direct numerical simulation of the relativistic 
magnetofluid. 

In the following report we will hrst discuss what is currently known about the 
dynamics of nonrelativistic MHD systems. We introduce the standard incompressible 
and compressible nonrelativistic MHD evolution equations as well as the ideal invariants 
for those systems. In section 4, we will introduce the relativistic MHD equations and the 
relativistic equivalents of the nonrelativistic ideal MHD invariants. We then describe our 
numerical experiment and present our results for a relativistic MHD code. We conclude 
by discussing the similarities and differences between the different systems. 

2. Nonrelativistic Incompressible MHD 

Work by Shebalin I25i m\, on ideal homogeneous incompressible MHD turbulence 
best demonstrates how the dynamics of a magnetofluid can differ from that of a 
hydrodynamic fluid. Plasma can be accurately modeled as a fluid made up of charged 
particles that are therefore affected by magnetic fields as well as particle-particle 
interactions. Because of this, the magnetic field becomes a dynamic variable in addition 
to density, pressure and the velocity of particles. For example. In MHD turbulence, an 
equipartition occurs and we expect kinetic and magnetic energy fluctuations to become 
roughly equal. Shebalin modeled the magnetofluid as a homogenous system where 
the same statistics are considered valid everywhere in the computational domain. He 
utilized periodic boundary conditions and spectral methods in order to study how the 
dynamics of different scales interacted without the addition of boundary errors. Much 
of his work focused on an ideal MHD system, where the magnetic and fluid dissipation 
terms were excluded. Below are the evolution equations used by Shebalin to describe 


the incompressible MHD system. 


o 
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Table 1. Invariants for Ideal Incompressible MHD. 


Case 

Mean Field 

Angular Velocity 

Invariants 

I 

0 

0 

Ej He, Hm 

II 

Bo ^0 

0 

E,He 

III 

0 

Do 7^ 0 

E, Hm 

IV 

Bq ^0 

Ho = 0 Bo 

E,Hp 

V 

Bo ^0 

Do 7 ^ 0 {Bo X Do 7 ^ 0 ) 

E 


By varying the mean magnetic field {Bq) and angnlar velocity (f2o) of the system, 
Shebalin was able to define five different cases with different invariants as shown in 
Table [H In snch a system there conld be as many as 3 ideal invariants; energy {E), 
and the psuedoscalars cross helicity {He) and magnetic helicity {Hm)- In addition, the 
invariant parallel helicity, Hp = He — ctHm {o' = VLq/Bq)^ can be formed from a linear 
combination of cross and magnetic helicity. In a hydrodynamic finid the ideal invariants 
are only energy and kinetic helicity. 

For an incompressible finid u(k,t) is the Fourier coefficient of turbulent velocity 
and b(k,t) is the Fourier coefficient of the turbulent magnetic field. The energy, cross 
helicity and magnetic helicity can be expressed in terms of these as: 

® =^E[i“Wt+iwr] (2“) 

k 

(2'>) 

k 

(2c) 

k 

A cubic computational domain with N grid points in each direction is assumed. The 
statistical mechanics of the system is defined by the Gaussian canonical probability 
density function (PDF): 

D = ^exp{-aE - /3He - '^Hm) 

Z = J exp{—aE — (3He — 'yHM)dT 

^ = [ ^DdT, ^ = ^ [ ^dt. 

Jr E Jq 

Where Z is the partition function and dF is the phase space volume. $ shows how to 
calculate the ensemble averages using the PDF while d> is the time average. If d) = $, 
the system is said to be ergodic but if d) 7 ^ $, it is non-ergodic. Here a, /3 and 7 are 
inverse temperatures. The ensemble average magnetic energy {Em) is always greater 
than or equal to ensemble average kinetic energy {Er), and the inverse temperature 


(за) 

( зб ) 

(3c) 
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terms can be found as a function of Em = 


2N^(f) 


(4 a) 


ct>{E - (t>) - HI 



0 

2(j)-E 


m 


7 =-^- a 


(4c) 


Phase portraits resulting from computer simulations of Shebalin's five cases show that 
coherent structures formed in many systems where the magnetofluid was experiencing 
turbulence [25l [26]. Coherent structure occurs when time-averaged physical variables 
in MHD turbulence have large mean values, rather than the zero mean values expected 
from theoretical ensemble predictions. MHD turbulence thus has broken ergodicity, 
which can be explained by hnding the eigenmodes of the system. One out of the 
four eigenvalues associated with each of the lowest wavenumbers will be very much 
smaller than the others; the eigenvariables associated with these very small eigenvalues 
grow to have very large energies compared to other eigenvariables; when this happens, 
an almost force-free state occurs in which large-energy eigenmodes are quasistationary 
while low-energy eigenmodes remain turbulent; thus, the predicted ergodicity has been 
dynamically broken. This is observed to occur even in dissipative systems because 
broken ergodicity in MHD turbulence manifests itself at the smallest wavenumbers 
(largest length scales) where dissipation is negligible, resulting in the ideal spectrum. 
In the case of ideal hydrodynamic turbulence, broken ergodicity can occur in a hnite 
model system, but only at the largest wavenumbers (smallest scales). When dissipation 
is added, the large wavenumber modes are most affected and their energy quickly decays 
away, so that broken ergodicity plays no role in decaying hydrodynamic turbulence. 

3. Nonrelativistic Compressible MHD 

Compressible MHD systems have not been studied as much as incompressible systems so 
here we will focus primarily on their invariants. We will assume that both incompressible 
and compressible systems share the same statistical mechanics and dynamics whenever 
the same invariants apply. In a nonrelativistic compressible MHD system; Energy and 
the Incompressible form of Cross Helicity are always conserved for a nondissipative 
system [H |22l [2^ (See Table 2). Compressible Cross Helicity, Hq = pHc, is not an 
ideal invariant in compressible MHD [T]. In the absence of a mean magnetic held 
and dissipation. Magnetic Helicity is also conserved [10]. The authors were unable to 
identify any literature showing the relationship between net angular velocity and ideal 
compressible MHD invariants. 

Given that our relativistic system is by default a compressible system, we naively 
expect to see that the same ideal invariants will apply for the relativistic system as 
the nonrelativistic compressible system. The equations for ideal compressible MHD are 
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Table 2. Invariants for Ideal Compressible MHD. 


Case 

Mean Field 

Invariants 

I 

0 

E, He, Hm 

H 

Po 7^ 0 

E,Hc 


similar to those of the incompressible system with the exception of the hrst equation. 

(5a) 

Vp- —B X CV X B) (5b) 

An 




dv 

p— = -p{v ■ Vv) 


OB 

dt 


= V X (n X B). 


(5c) 


4. Relativistic MHD Systems 


The fluid and electromagnetic components of the relativistic MHD equations are 
developed from several well-known equations [ 8 ]. They include the conservation of 
particle number, the continuity equation, the conservation of energy-momentum, the 
magnetic constraint equation and the magnetic induction equation. For a system 
consisting of a perfect fluid and an electromagnetic held, the ideal MHD stress-energy 
tensor is given by 


b 2 

= (poh + + (P + 

P = 6po(f - 1) 

P 

h = 1 -b e H- 

Po 

b^ = -^Bt 


\/4^ 


(n) 


BU = 


( 6 a) 

(65) 

(6c) 

(6<i) 

(6e) 


BU = \(- + Blp 


( 6 /) 


Here, P is the huid pressure, po is density, P* is magnetic held, is four-velocity, h 
is the enthalpy, e is specihc internal energy, and 6 ^ is the magnitude of the magnetic 
vector held squared. We dehne pressure in terms of the energy density using the T law 
equation of state with, T = |. P = (T — l)pe at most energies and P oc p'" at extremely 
low energies. The evolution equations where given by Duez as [5]: 


dtp* = -dj{p*v^) 

dtf = - p*P) -b s 


(7a) 

(76) 
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(7c) 

- PB^). 

(7<i) 


Here, p* is conserved mass density, f relates to energy density. Si is momentum density, 
is related to the magnetic field and s is the source term, 'jij is the three metric and 
a is a lapse term related to the time evolution of the simulation. The determinate of 
the three metric and lapse are both set to unity because we are using the Minkowski 
metric and Geodesic slicing conditions. 



= a^pou^ 

(So) 

f 

= - p. 

m 


= = [p^h + au^y/^l?)ui - 

( 8 c) 



( 8 d) 


Notice that unlike the nonrelativistic system, we use the stress-energy tensor within the 
evolution equations so that 4-momentum conservation is built into the system. This 
results in a set of equations that look very different from that of the nonrelativistic 
system. According to work by Yoshida [2H1EO] , in addition to 4-momentum, relativistic 
systems are expected to conserve a quantity called Relativistic Helicity. It is defined 
below using the canonical 4-momentum density, V, of the system. 

= (p . (V X V),Vo{V xV)+Vx (VPo + W) (9) 

Here the canonical 4-momentum density V = {Vo,V) is a combination of mechanical 
and electromagnetic momentum densities. The conservation of Relativistic Helicity is 
then effectively, J = 0. The canonical 4-momentum can be expressed as the 

sum of mechanical and electromagnetic momentum, If we ignore the 

electromagnetic fields, we recover a relativistic version of Cross Helicity Density, kc- If 
we set the particle’s mechanical momentum to zero, we recover a relativistic version of 
Magnetic Helicity Density, km- 

5. Methodology 

For this numerical experiment, we calculate Energy Density (E), Relativistic Helicity 
Density (k). Cross Helicity Density (kc); and Magnetic Helicity Density (km), 
numerically using the code described later in this section. These variables are defined 


as shown below. 

B (10=) 

Po = p*huo - P (106) 

Pi = p^hui ( 10 c) 

n =(S-(VxS),(r + p,)(VxS) + Sx (V(r + p*) + doS)) (lOd) 
Kc = (P ■ B, PqB — P X T) (lOe) 

KM = iA-B,AoB~AxS) ( 10 /) 
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Here the magnetic field is related to the vector potential by the equation, B = V x A. 
The electric field is defined using the MHD conditions, 8 = B x w. We can test to 
see if the Helicities are invariant by comparing numerically calculated time derivatives 
of lA to their predicted value at each time-step using the equation OqiA + din} = 0 . 
We normalize the result using the L2 norm of the calculated divergence so all results 
are on the same relative scale. We then integrate the result over the volume of our 
computational domain. If the normalized error is dominated by the truncation and 
round-off errors, we can assume that the system is invariant. For the normalized error 
in energy we simply look at the difference in energy at two different time levels divided 
by the total energy at that time level. 


Error K = 


Error p, = 


E{t) - E{t - At) 

W) 


d^x 


(11a) 

(lib) 


In order to study the invariants of the relativistic MHD equations, we used a code 
called FixedCosmo which was originally written by one of us |9] to study the dynamics 
of primordial plasma turbulence. This code was developed using the open source Cactus 
framework (www.cactuscode.org). Cactus was originally developed to perform numerical 
relativistic simulations of colliding black holes but it’s modular design has since allowed 
it to be used for a variety of Physics, Engineering and Computer Science applications. 
It is currently being maintained by the Center for Computation and Technology at 
Louisiana Sate University. Cactus codes are composed of a flesh (which provides the 
framework) and the thorns (which provide the physics). FixedCosmo is a collection of 
thorns. It uses the form of the Relativistic MHD equations described by Duez jH] and is 
written in a combination of F77, F90, C and C-I--I-. This code is parallelized and capable 
of using several different differencing methods such as second order finite differencing, 
fourth order finite differencing and spectral differencing. Although the code is capable 
of utilizing artificial viscosity and HRSC, neither was used for this project. 

Because the objective of this study is to test the ideal relativistic MHD system, we 
complete a series of runs in a “high-energy” regime. The parameters used approximate 
that of the early universe around the electroweak scale. This is done so we can apply 
the results to any relativistic MHD system. Table 3 shows a matrix of the test runs. 

Each data run utilized Fourier spectral differencing on a grid with 64 x 64 x 64 
internal data points. We ran these simulations for about 7500 iterations or over 10“® 
s of physical time. The electron oscillation time for the “high energy” regime is about 
1.1 X 10“^^ s so the simulations appear to have run long enough to witness the full 
dynamics of the system. The simulation domains where all set to 4m x 4m x 4m. Also, 
the code utilizes geometerized units {c = G = h = 1) so parameters where translated 
from SI units to units of length for the use of the calculations and then back to SI units 
for the output. Time is therefore translated into seconds by dividing the output time 
by the speed of light. 
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Table 3. High Energy Numerical Simulations 


Variables 

Gase 1 

Gase 2 

Gase 3 

Gase 4 

Gase 5 

Max Velocity (c) 

0.25 

0.25 

0.25 

0.25 

0.25 

hlj; (rad/s) 

0 

0 

0 

0 

0 

VLy (rad/s) 

0 

0 

0 

0 

0 

(rad/s) 

0 

0 

0.35 

0.35 

0.35 

Init Temperature (K) 

2.8el5 

2.8el5 

2.8el5 

2.8el5 

2.8el5 

Init Density {kg/m^) 

9.7e29 

9.7e29 

9.7e29 

9.7e29 

9.7e29 

Max Magnetic Field (G) 

1.0el3 

1.0el3 

1.0el3 

1.0el3 

1.0el3 

(G) 

0 

2.0el3 

0 

0 

2.0el3 

By (G) 

0 

0 

0 

0 

0 

B, (G) 

0 

0 

0 

2.0el3 

0 


6. Results 

Truncation errors were found by doubling the resolution and measuring the change in 
the observed total errors. By assuming that the Euler Method, used to calculate 
is hrst order and that the numerical errors and variations from exactly conserved values 
are additive, the truncation errors can be calculated from 

\Etruncatioin\ (12) 

If the truncation errors are within an order of magnitude of the normalized errors, we 
can conclude that the system is invariant. We also found it impossible to completely 
eliminate the mean magnetic held and mean angular momentum in all cases. A mean 
magnetic held (on the order of 1% of the maximum held) remained in every case. Also, 
each case seemed to have a small angular velocity, also less than 1% of the huid velocities 
within the simulation. The authors feel that these residual quantities where not enough 
to signihcantly disrupt the system and could be safely ignored. 


Table 4. High Energy Simulation Results 


Mean Magnitude of Errors 

E 

K 

i^c 

Km 

Gase 1 

1.8e-9 

2.8e-14 

9.6e-5 

l.Oe-2 

Gase 2 

1.8e-9 

2.8e-14 

9.7e-5 

2 .8e-2 

Gase 3 

8 .1e-10 

0.0 

6.6e-5 

7.1e-3 

Gase 4 

7.7e-10 

0.0 

6.0e-5 

1.9e-2 

Gase 5 

7.8e-10 

0.0 

6.6e-5 

1 .8e-2 

Truncation Errors 

2.1e-9 

2.2e-14 

2.9e-5 

2.0e-3 


The results in Table 4 were calculated by averaging the absolute values of 
normalized errors. As one can see. Energy and Relativistic Helicity appear to be 
conserved in every case given that the normalized error appears to be dominated by 
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truncation errors for both. Cross Helicity may be conserved in every case since the 
calculated truncation error appears to be within an order of normalized errors in every 
case. Magnetic Helicity does not appear to be conserved in any of the cases. Normalized 
Errors for Cross Helicity Conservation appear smaller in cases where a large mean 
angular velocity is present. Deviations in Magnetic Helicity Conservation are smallest 
in the absence of a large mean magnetic held. 

7. Discussion 

Our results show that in the high-energy Relativistic MHD regime only Energy and 
Relativistic Helicity are clearly conserved. We are not able to conclusively prove Cross 
Helicity conservation. Magnetic Helicity conservation is questionable in this system. 
This is not an unexpected result but it does raise several interesting questions which 
lie beyond the scope of this article. Does the potential lack of Cross and Magnetic 
Helicity Conservation effect the dynamics of the relativistic system when it comes to 
phenomena such as inverse Energy Cascade or the Kolmogorov Energy Spectrum? How 
do magnetic dynamos in relativistic MHD systems function? Are there any other 
overlooked dynamics in relativistic MHD systems? These are all questions which we 
hope to address in future numerical studies. 
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